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We present the results of our Monte Carlo simulation of the Ising-0(3) model on two-dimensional (2D) and 
quasi-2D lattices. This model is an effective classical model for the stacked square-lattice J1-J2 Heisenberg 
model, where the nearest-neighbor (Ji) and next-nearest-neighbor (J 2 ) couplings are frustrated and we assume 
that Ji is dominant. We find an Ising ordered phase in which the 0(3) spins remain disordered in a moderate 
quasi-2D region. There is a single first-order transition for a sufficiently large 3D coupling, in agreement with 
a renormalization group treatment. The subtle region in which the single transition splits into two transitions 
is also discussed and compared against recent measurements of two very close transitions in BaFe2As2. Our 
results can provide a qualitative explanation of the experiments on ferropnictides, namely the observed sequence 
and orders of the structural and magnetic transitions, in terms of the ratio between the inter-layer and intra-layer 
coupling. 
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I. INTRODUCTION 

The dimensional crossover [from one-dimensional (ID) or 
2D phenomena to 3D ones] has been an important issue in 
condensed matter physics because there are several quasi-low- 
dimensional materials that exhibit complex behaviors. En- 
hanced quantum fluctuations in low-dimensional systems can 
give rise to novel quantum states of matter. Although a small 
3D coupling always exists in real systems and ordered phases 
are usually stabilized at low temperatures, the dimensional 
crossover is still governed by the low-dimensional physics. 

When a system combines fluctuating continuous and dis- 
crete degrees of freedom, the interplay between them leads 
to an even richer dimensional crossover. The simple rea- 
son is that they respond to a weak 3D coupling in a quali- 
tatively different way. This physics has recently attracted par- 
ticular interest in the context of parent compounds of iron- 
based superconductors.- - - Among the iron-based compounds, 
the 1111-type quasi-2D materials i^FeAsO (R is a rare earth 
ion) 4,5 and the 122-type 3D compounds AFe2As2 (A is an 
alkaline earth ion) 6-9 constitute a subgroup with the follow- 
ing low-energy properties: 1 (i) they are metallic, (ii) they un- 
dergo a tetragonal-orthorhombic structural transition, and (iii) 
they stabilize a stripe-like spin-density-wave (SDW) order at 
a wavevector (7r, 0) at the lowest temperatures. The struc- 
tural transition is an Ising-like transition in the sense that it 
breaks a Z2 spatial symmetry. The continuous and discrete 
characters of these two broken symmetries lead to qualita- 
tively different magnetic and structural fluctuations whose in- 
terplay should be strongly affected by the magnitude of the 
3D coupling. Indeed, experiments have confirmed that the 
structural and magnetic transitions take place simultaneously 
via a first-order transition in most of the nearly 3D undoped 
122 compounds whereas the lattice distortion occurs at a 
slightly higher temperature in the more quasi-2D 1111 mate- 
rials (Fig. [T]). 5 In the latter case, both transitions seem to be 
of second order or at least very weakly first order.- The prox- 
imity between the two transitions suggests that the magnetic 
ordering plays a central role in the lattice distortion*!^ Previ- 



ous density functional studies also indicate that the structural 
distortion may be driven by the interaction between magnetic 
degrees of freedom^ 

These qualitative and ubiquitous features of the iron pnic- 
tides are expected to be universal in the sense that they should 
only depend on symmetry, dimensionality, number of com- 
ponents of the order parameter, and range of interactions. 
For this reason, it has been argued that the J\- J 2 Heisenberg 
model 13 is very useful in spite of its simplicity as a purely lo- 
cal spin modeL 10,11 Here, we consider this model on a stacked 
square lattice. The model includes in-plane nearest-neighbor 
(Ji) and next-nearest-neighbor {J 2) exchange couplings that 
lead to a high degree of frustration. In addition, we assume 
a weak unfrustrated inter-layer coupling. It is known that 
this model leads to a stripe-like antiferromagnetic ordering 
with wavevectors (ir, 0, q z ) or (0, 7r, q z ) for sufficiently large 
J2/J1' Here q z = 0, tt depends on the sign of the inter-layer 
coupling. In d = 2, this collinear order becomes stable for 
J2/J1 > 0.66 in the quantum limit (S = 1/2) 14 and for 
J2/J1 > 0-5 in the classical limit (S — » 00). The J1-J2 
Heisenberg model can be derived from the multiband Hub- 
bard model proposed for the iron-based compounds by taking 
the strong-coupling limit— 

A key observation here is that the stripe-like magnetic or- 
dering breaks a discrete Z2 symmetry, associated with two 
possible bond orderings, in addition to the continuous SU(2) 
symmetry. Therefore, in addition to the spin- wave excitations 
associated with the broken continuous symmetry, the model 
with dominant J2 also includes low-energy Ising-like degrees 
of freedom as pointed out by Chandra et al. 16 Here, the dis- 
crete Z2 symmetry corresponds to a tt/2 rotation of the square 
lattice. The bond ordering of the Ising-like degrees of free- 
dom triggers the tetragonal-orthorhombic structural distortion 
in the presence of a finite spin-lattice coupling. 17 For this rea- 
son, the Ising-like ordering is usually interpreted as a struc- 
tural transition when referring to real compounds J^ii We will 
denote the transition temperatures for the structural distortion 
and the magnetic ordering by T c \ and T C 2, respectively. Since 
the magnetic ordering cannot exist without the Ising order- 
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FIG. 1. (Color online) Schematic pictures of the sequence of the 
transition(s) observed in the parent compounds of the iron-based su- 
perconductors. Shown are the disordered state with tetragonal lattice 
symmetry at high T's (left), the Ising ordered state with broken lat- 
tice symmetry and short-range magnetic ordering at intermediate T 
(center), and the lowest-T state with broken lattice symmetry and 
long-range stripe-like magnetic ordering (right). The symmetry also 
allows for the other lattice distortion obtained by rotating the pre- 
sented lattice by tt/2. 



ing, it is clear that T& > T C 2 10,11 in agreement with exper- 
imental observations. It is then natural to ask what are the 
additional experimental aspects that can be explained quali- 
tatively with a local spin model. Finally, the study that we 
present in this paper is even more relevant for insulating com- 
pounds such as Li2VOSi04 and Li2VOGe04^ These com- 
pounds are believed to be well described by the J1-J2 model, 
and an experimental signature of the Ising-like structural tran- 
sition just above the magnetic transition has been suggested 
for Li 2 VOSi0 4 ^ 

Previous theoretical treatments of the dimensional 
crossover in the J1-J2 Heisenberg model at T > relied 
on approximate methods such as large- N expansions^ a 
random-phase approximation (RPA) (or a layer mean-field 
theory)ri or a phenomenological Landau mean-field theory. 19 
Although all of these treatments agree on that the model has 
an Ising ordered phase with unbroken spin symmetry in a cer- 
tain quasi-2D region, the precise form of the phase diagram 
as a function of the inter-layer coupling is still unknown. 
Moreover, the results obtained by these approximations have 
several contradicting points relative to the detailed structure 
of the phase diagram, especially when the inter-layer coupling 
becomes stronger and the two transition temperatures become 
closer to each other. The main motivation of the present study 
is to resolve these contradictions by applying a controlled 
numerical method. We introduce an unfrustrated classical 
effective model suitable for the J1-J2 Heisenberg model 
with dominant J 2 • By using a classical Monte Carlo method, 
we show that the interplay between the Ising and magnetic 
degrees of freedom leads to a first-order transition when the 
magnitude of the inter-layer coupling is sufficiently large. 
We also provide a renormalization group (RG) argument 
supporting this observation. Finally, we numerically identify 
the Ising ordered phase in the quasi-2D region and present the 
corresponding phase diagram. The subtle region where the 
single transition splits into two transitions is also discussed 



and compared against recent measurements of two very close 
transitions in BaFe2As2*22i2i 



II. MODEL 

We start by considering an unfrustrated classical model that 
is an effective Hamiltonian for describing the physics of the 
J 1- J 2 Heisenberg model with dominant J 2 on a quasi-2D sys- 
tem of weakly coupled square-lattice layers. The Hamiltonian 
associated with this so-called Ising-0(3) modelii^ is 



H 



(id) 



Jij (1 + CFiVj) Si • Sj 



(1) 



where G{ and denote the classical Ising and 0(3) spins, 
respectively. The spatially anisotropic coupling constant is: 

j I J if (i, j) are on the same layer, 
lJ 1 J z if (i, j) are on the nearest neighbor layers, 

where < J z < J. In the following we use J = 1 as the unit 
of energy unless otherwise specified. 

As we mentioned above, when J2 > Ji/2, the J1-J2 
Heisenberg model shows magnetic orderings with wavevec- 
tors (7r, 0, q z ) or (0, tt, q z ). The stripe-like ordering breaks the 
lattice rotational symmetry as well as the 0(3) spin symmetry. 
The ordered state consists of two interpenetrating a/2 x a/2 
sublattices, each of which shows a simple Neel order, and 
the inter- sublattice coupling exactly vanishes in the absence 
of thermal and quantum fluctuations. Fluctuations stabilize 
the stripe-like order via the generation of a biquadratic cou- 
pling between the two sublattice order parameters favoring the 
collinear (stripe-like) spin configuration. 16 This is a clear ex- 
ample of order by disorder*i&22r— The effective Hamiltonian 
that describes the transition to this broken symmetry state in 
the long- wavelength limit is: 11 

H GS = fd d x f^|V0a| 2 + r|0 a | 2 +ii|0 a | 

J [ a—A,B ^ Z 



■ u A b\4>a\ 2 \4>b\ 2 + A (4>a • 0b) 



• (3) 



Here, 4>a and 0# represent the three-component sublattice 
magnetic order parameters. The first three terms describe 
intra- sublattice fluctuations, while the other two quartic terms 
describe the inter- sublattice couplings allowed by symmetry. 
u and uab are positive, whereas A is negative. The negative 
A forces 4>a and 4>b to be collinear and {4>a • 4>b) becomes 
the Ising-like order parameter that decides whether the stripe- 
like spin configuration is "vertical", as in Fig.Q] or "horizon- 
tal." The model can in principle retain the Ising ordered phase 
while the magnetic ordering is only of short range. In this 
phase, we have (4>a) = (4>b) = but (4>a • 4>b) ^ 0. 16 

The Ising-0(3) Hamiltonian is an effective model for the 
J1-J2 Heisenberg model in the sense that both Hamiltoni- 
ans are described by the same effective theory in the long- 
wavelength limit if J2/J1 is sufficiently large. The corre- 
sponding derivation requires introducing auxiliary 0(3) fields 
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4>a ^ S and 4>b ^ crS, as is described in Appendix lAl The 
Ising variables, (c^), of the Ising-0(3) model correspond to 
((f) a • 4>b), while (Si) describes the local magnetic ordering. 
More heuristically, we can first write down a two-sublattice 
classical spin model of the form 27 

^coupled = — Jij (^A,i ' &A,j + &B,i ' &B,j) 

+ ^A(S^.Sb,,) 2 , (4) 

and take the strong biquadratic coupling limit A — >> — oo. 
Based on symmetry arguments, it is clear that the long- 
wavelength limit of i^coupied is also described by H Q ff. In the 
limit A — > — ex), one can write S^,i = Si and Sb,i = 
and this procedure yields the Ising-0(3) model. The Ising-XY 
model, in which S represents a classical XY spin, was stud- 
ied extensively as an effective model of the fully frustrated 
XY Hamiltonianr&~— and an arbitrary coefficient was some- 
times included in front of a^aj in Eq. (Q]). This prefactor must 
be equal to unity in the present case because the original J\- 
J2 model is invariant under exchange of the two sublattices 
A ±+B. 

The factor (1 + (Ji(Jj) in Eq. (Q} can be viewed as the 
effective coupling of the nearest-neighbor classical 0(3) spins 
Si and Sj . Because it is non-negative, ferromagnetic align- 
ment of the 0(3) spins is always favored. This in turn implies 
that the mean-field coupling, (Si • Sj), between the Ising 
variables is also ferromagnetic. The form of (1 + cr^-) 
also implies the following restriction on the ordering of the 
0(3) spins. Let us consider a case in which the Ising spins are 
disordered and we divide the system into ferromagnetic clus- 
ters of the Ising variables. We will now quench the Ising spin 
configuration and perform a partial trace on the 0(3) spins. 
We see that an 0(3) spin of a given cluster cannot correlate 
with 0(3) spins in different clusters because (1 + Oi<Jj) 
vanishes at the boundary between the clusters, i.e., on bonds 
between anti-parallel Ising spins. Consequently, the 0(3) 
spins cannot order without Ising ordering, because a perco- 
lating Ising cluster is required to have 0(3) ordering. In this 
way, we arrive again at the general inequality T c \ > T C 2. 

In the 2D limit, J z = 0, we expect a finite-temperature 
Ising transition, while the 0(3) spins must remain disordered 
at any finite temperature because of their non-invariance un- 
der a continuous symmetry. 31 The Ising transition was con- 
firmed in the 2D classical J1-J2 Heisenberg model on the 
square lattice by a Monte Carlo simulation^ In d = 2, we 
expect that the behavior of the 0(3) spins far below the Ising 
transition should be very close to that of the square-lattice 
Heisenberg model. A finite value of J z induces a dimen- 
sional crossover: the Ising transition is shifted to a higher 
temperature, and the 0(3) spins also become ordered at a low 
enough value of T. A simple RPA argument predicts that 
T c i — ~ (Jz) 1 ^, where 7 is the 2D Ising exponent 
and T C 2 ~ — 47r/ In J Z M^~— The qualitative difference orig- 
inates in the different critical behaviors in the 2D limit. The 
power-law correlations of the Ising variables at T = T c \ are 
qualitatively different from the essential singularity at T = 



for the 2D 0(3) model. 

What should we expect well inside the 3D regime that oc- 
curs for large enough J z l The RPA argument, which de- 
scribes different order parameters independently, is inappro- 
priate in this regime because the interplay of the Ising and 
0(3) variables is expected to be strong. A previous large- A 7 ' 
treatment indicates that the two transitions never merge and 
remain of second order A phenomenological Landau mean- 
field theory predicts that the 0(3) transition should become 
of first order before merging with the second-order Ising tran- 
sition. 19 In addition, this theory predicts a single first-order 
transition for intermediate values of J z and a single second- 
order transition for larger values. 

On the other hand, an RG analysis suggests that the merged 
transition will always be of first order. The reason is that a 
one-loop epsilon expansion applied to a generalization of H t a 
in a different context (amorphous magnets) leads to no sta- 
ble fixed points, 36 implying the absence of the scale invari- 
ance that is characteristic of second-order transitions (see Ap- 
pendix [B] for details). This RG result contradicts the above- 
mentioned Landau theory that allows for a single-second or- 
der phase transition when the inter-layer coupling J z is large 
enough^ 



III. METHOD 

We simulated the Ising-0(3) model on the square lattice 
(J z = 0) and the quasi-2D anisotropic cubic lattice (0 < 
J z < 1) using the Monte Carlo method. We employ a cluster 
Monte Carlo method in which the clusters of Ising and 0(3) 
spins are updated alternatively. Updates of the 0(3) spins take 
place for fixed Ising variables based on the Wolff algorithm 37 
and the quenched coupling constant (1 + <Ji<Jj). The same 
idea is used for updating the Ising spins, with JijSi ■ Sj play- 
ing the role of the effective coupling. The main difference 
is that J^Si • Sj can be negative even though its average is 
positive. This fluctuation in the sign of the effective coupling 
can dynamically induce an effective frustration that makes the 
cluster updates less efficient. In practice, it turns out that 
the fluctuating sign effect does not matter for nearly spatially 
isotropic systems in d = 2 and 3, but it does matter when J z 
is small. To remedy this problem, we incorporate the clus- 
ter update with local and semi-global updates of Ising spins. 
In semi-global updates, we allow the clusters to expand only 
within a given layer by employing the Swendsen- Wang-type 
multicluster scheme. The Boltzmann weight factor related to 
the inter-layer couplings is absorbed in the cluster-flip proba- 
bility to satisfy detailed balance. Relatively small intra-layer 
clusters are expected to be flipped in a collective way by the 
application of this empirical method. 



IV. RESULTS 

We now show the results of our Monte Carlo simulation. 
The Ising-0(3) model exhibits (i) a second-order transition in 
the Ising universality class for J z = 0, (ii) a single first-order 
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FIG. 2. (Color online) Phase diagram of the quasi-2D Ising-0(3) 
model. The Ising ordered phase ((S) = 0, (a) / 0) appears in the 
quasi-2D region J z < 0.01. For large enough inter-layer coupling, a 
single first-order transition separates the paramagnetic state ((S) = 
0, (a) = 0) from the lowest-T phase that simultaneously breaks the 
Z2 and 0(3) symmetries ((S) 7^ 0, (a) 7^ 0). The detailed structure 
for 0.01 < J z < 0.0204 remains to be clarified; see the text. The 
line representing a phase boundary between the lowest-T phase and 
the Ising ordered phase is a schematic one. The other lines are guides 
to the eye. 



transition when J z is sufficiently large, and (iii) two second- 
order transitions in a moderate quasi-2D region ( J z < 0.01). 
The corresponding phase diagram is presented in Fig. O We 
emphasize that the two transitions merge into a single first- 
order transition that simultaneously breaks the Z 2 and 0(3) 
symmetries, in agreement with the RG treatment. As for the 
split transitions in the quasi-2D region, our numerical results 
suggest that the transitions are in the 3D Ising and the 3D 0(3) 
universality classes. There is still some level of uncertainty 
relative to the merging of the Ising and 0(3) transitions around 
0.01 < J z < 0.0204. However, our results show some subtle 
features that agree with the Landau mean-field theory in this 
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region^ 



A. Ising transition of the 2D model 

We begin with the 2D Ising-0(3) model. The maximum 
lattice size that we studied is L = 192. Due to the effective 
ferromagnetic coupling between Ising variables generated by 
the nearest-neighbor 0(3) spins, the 2D system undergoes a 
finite-T transition where only the Ising variables are critical 
while the 0(3) spins remain disordered. This phase transition 
can be analyzed very efficiently by introducing the following 
dimensionless scaling parameters: the Ising Binder parameter 
defined by U a = (<T 4 )/(a 2 ) 2 , with a = L~ 2 a u and the 
second moment correlation length 39 of the Ising variables in 
units of L, £cr/£, with 
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FIG. 3. (Color online) (a) £ a /L and (b) U a of the 2D Ising- 
0(3) model. The lines are guides to the eye. The horizontal lines 
marked as "2D Ising (CFT)" indicate the universal critical values of 
the square-lattice Ising model given in Ref . 38|. The inset shows their 
FSS plots for L > 48, where v = 1 and T C T = 1.0514(3). 



Here <r q is the Fourier mode at the lowest nonzero momentum, 
q = (27r/L, 0) or (0, 2tt/L), for a given lattice. These vari- 
ables asymptotically cross each other at the Ising transition 
point as a function of the system size. As shown in Fig. [3j 
the crossing point found for the larger lattices leads to a value 
of TlY = 1.0514(3). It is also known that the values of these 
scaling parameters should be universal at their crossing points. 



CI 

a 




1/v {rv ^Dx , X 2D 
L U" A C 1 ) 1 A cl 



4. (Color online) FSS of the Ising correlation function 
(l/2,l/2) at m e largest distance in a given 2D system defined by 
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FIG. 5. (Color online) The inverse susceptibility of the 0(3) order 
parameter of the 2D Ising-0(3) model (J z = 0). The vertical line 
indicates the 2D Ising transition temperature T^f . 



In agreement with this expectation, our results are consistent 
with the universal values for the 2D Ising model, which are 
obtained by using the conformal field theory (see Fig. [3]) 3 
Our final confirmation of the 2D Ising universality class is a 
finite-size scaling (FSS) plot of the dimensionless quantities 
shown in the inset of Fig. [3j and another FSS plot of the cor- 
relation function for the Ising spins at the most distant sites of 
a given system size (see Fig. [4]), 

G(l/2,l/2) = (Wi) with r y = (L/2, L/2) . (6) 
We assume 77 = 0.25 and v = 1 in these FSS plots. 

B. Separate transitions in the quasi-2D regime 

We now discuss the quasi-2D regime. The finite inter-layer 
coupling makes the Ising order three-dimensional and it also 
stabilizes the 0(3) spin order at T > 0. Of our particu- 
lar interest, motivated by the separate transitions in the 1111 
compounds, is the region where the Ising and the 0(3) transi- 
tions occur at different temperatures. The order of magnitude 
of the inter-layer coupling in such a region can be estimated 
from a simple RPA argument ,1122^2 The previous RPA esti- 
mates of T c i and T C 2 were obtained from the following con- 
ditions: J zX f{T cl ) ~ 1 and J zX ™{T c2 ) ~ 1 where X f 
and Xm are m e susceptibilities of the Ising and 0(3) variables 
for J z =0. To obtain a rough estimate of J z in the regime 
of separate transitions, we define J* as the magnitude of the 
inter-layer coupling for which the RPA estimate of T C 2 coin- 



cides with T c 2 f 



1.0514(3), i.e., (J*)" 



2D 



m- By 



using our numerical estimation of Xm shown in Fig.[5j we ob- 
tain J* = 0.00765(2). Figured in conjunction with the RPA 
argument, also suggests that J z < 0.001 is required to obtain 
a separation T c \ — T C 2 of order 0.01 J. 

Unfortunately, a direct finite-size study for J z < 0.001 is 
not simple because the highly anisotropic correlations would 
require very large system sizes even for simulating a few lay- 
ers (see, e.g., Ref. 35). For this reason, we concentrate on 
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FIG. 6. (Color online) Two separate transitions for J z — 0.01: 
(a) specific heat, (b) the Binder parameters, (c) the normalized cor- 
relation lengths in the x direction, and (d) the normalized correlation 
lengths in the z direction. The vertical lines show the estimated crit- 
ical temperatures based on the FSS analysis. The lines are guides to 
the eye. 



J z = 0.01. (We also simulated the model for J z = 0.006667 
and observed essentially the same phenomena.) Although 
this value of J z is slightly larger than J*, the two transi- 
tions still occur at different temperatures. To reduce the finite- 
size effects induced by the spatial anisotropy, the system is 
set to have a tetragonal shape L x x L y x L z with periodic 
boundary conditions (L x = L y = L). The aspect ratio 
r = L z /L = 1/12 was determined in such a way that 
(La/2, 0,0) « (0,0,2^/2) and G m (Z^/2, 0, 0) w 
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FIG. 7. (Color online) FSS plots for the higher-T Ising transition 
for J z = 0.01 [ti = (T - Tci) /Tci]. The upper panel (a) shows 
the FSS of U a . The inset shows the relation to the real temperature 
T; the thick line indicates the lower-T transition close to which a 
severe finite- size effect appears. The middle panel (b) shows the FSS 
of £%/L x , with the inset showing the FSS of £*/L z . The lower 
panel (c) shows the FSS of G a Lx/2 = G a (L x /2, 0, 0) and G a Lz/2 = 
G a (0, 0, L z /2) (the inset) [see Eq. Q]. 



FIG. 8. (Color online) FSS plots for the lower-T 0(3) transition 
for J z = 0.01 [t 2 = (T — T c2 ) /T c2 \. The upper panel (a) shows 
the FSS of Um- The inset shows the relation to the real temperature 
T; the thick line indicates the higher-T transition close to which a 
severe finite- size effect appears. The middle panel (b) shows the FSS 
of ^n/Lx, with the inset showing the FSS of ^/L z . The lower 
panel (c) shows the FSS of G^ x/2 = G m (L x /2, 0, 0) and G^ z/2 = 
G m (0, 0, L z /2) (the inset) [see Eq. Q]. 



G m (0,0,T 2 /2), where 

G a (r x ,r y ,r z ) = (a l a J ), G m (r x , r y , r z ) = (S< • S,) (7) 

with Yij (r r z ) are the correlation functions of the 
Ising and the 0(3) spins. We studied a range of system sizes 
from L = 18 to L = 144. 

In Fig. [(Ja), we show the specific heat per site, C = 
N~ 1 f3 2 (( K H 2 ) - (H) 2 ), where N = rL 3 is the number of 
sites. The double-peak structure of the specific heat with in- 
creasing depth for larger systems gives a first indication of 
two separate transitions. This observation is supported by 
the behavior of the dimensionless scaling parameters. Here 
we use the Binder parameters, U a = (a 4 ) /(a 2 ) 2 with a = 



N' 1 ^ and U m = (m 4 )/(m 2 ) 2 with m = N' 1 £\ S,, 
and the normalized correlation lengths, / T M and of 
the Ising and the 0(3) order parameters along the intra-layer 
(/jl = x) and the inter-layer (fi = z) directions. ^ and £ m are 
defined by: 



^ = 



(* 2 )/(<)- 

\ 4sin 2 (7r/L M ) 



\ 



n| 2 V(|m q J 2 



4 sin 2 (n/Lf,) 



(8) 



(9) 



where q M is the lowest nonzero momentum for a given lattice 
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in the /i direction. As shown in Figs.[6lb)-[6jd), the curves of 
these quantities exhibit a clear tendency to intersect for larger 
lattices, and the crossing temperature varies depending on the 
order parameter to which they are related. While we find no 
peak in U a , there is a peak structure in U m for small systems, 
but this is suppressed for larger lattices. The non-divergent be- 
havior of the Binder parameters eliminates the possibility of 
strongly first-order transitions.— The FSS plots are presented 
in Figs. [7] and [8] Based on symmetry arguments, we assume 
that the higher-T (lower-T) transition is in the 3D Ising [0(3)] 
universality class, and we use the corresponding critical expo- 
nents that are available in the literature: r]i s = 0.03639(15) 
and V[ s = 0.63012(16) for the 3D Ising universality class^- 
and 7? H = 0.0375(5) and m = 0.7112(5) for the 3D 0(3) 
universality class.— Although finite-size effects are still se- 
vere for the explored system sizes, we can see an asymptotic 
tendency toward data collapse. This observation supports the 
assumed universality classes, which leads to T c \ = 1.0610(7) 
and T C 2 = 1.0575(2). The rather small separation is indeed 
expected because J z = 0.01 is relatively large in comparison 
with J*. 

As for the sizable sub-leading corrections observed in these 
scaling plots, it appears that they are largely due to the prox- 
imity of the two transitions and/or the spatial anisotropy of 
the correlations. Naturally, the proximity effect is expected 
to appear in the low- (high-) temperature side of the scaling 
plots for the Ising [0(3)] transition. For instance, a devia- 
tion appears in the scaling plots around the 0(3) transition for 
l i/vk ( T _ Tc2 ) / Tc2 > i an d L = 72 (Fig.©. This devi- 
ation is most likely due to the proximity to the Ising transi- 
tion. A similar effect is observed in the FSS plots of the Ising- 
like transition for L x l Vls (T - T cl ) /T cl < -2 and L = 72 
(Fig. [T]). However, these finite- size effects disappear rapidly 
for larger system-sizes. As for the finite-size effect due to 
spatial anisotropy, we find that the crossing value of / 
depends on /i = x, z, although the crossing value of ^ m /L^ is 
almost independent of fi [Figs.^c) and[6td)]. This observa- 
tion suggests that the aspect ratio r = 1/12 for J z = 0.01 
is appropriately tuned to investigate the 0(3) transition of 
this system, but it is not perfectly tuned for investigating the 
Ising transition. The observation of £ a /L x > / L z implies 
a shortness of the effective inter-layer coupling of the Ising 
spins in the simulated finite system (with r = 1/12) and such 
an anisotropy effect might produce sub-leading corrections to 
the scaling behavior near T = T c ±. However, we believe that 
these corrections will not affect our conclusions significantly. 
In particular, our conclusion about the separation of the tran- 
sitions for J z = 0.01 does not change. The reason is that our 
estimation of T C 2 is accurate enough and a modified aspect ra- 
tio of r < 1/12 will never lower the estimation of T c \ because 
it has the effect of enhancing the effective inter-layer coupling 
between Ising variables. 



C. 3D system with a large inter-layer coupling 

We now discuss the 3D regime, where the inter-layer cou- 
pling is sufficiently large. We note that previous analytical 
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J Z =0.0625,T=1.1232,L=24 ^ 
J Z =0.04,T=1.0969,L=40 i-= 
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FIG. 9. (Color online) Bimodal internal energy density distribution 
at the first order phase transition. The inset shows the J z -dependence 
of the peak-to-peak distance AE of the distributions. 
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(c) J z = 0.0204 
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FIG. 10. (Color online) Binder parameters U m and U a near the 
first order transition for (a) J z = 0.0625, (b) J z = 0.0278 and (c) 
J z = 0.0204. 



studies suggested different scenarios in this region. A large- 
N approximation 10 predicted a scenario with two separate 
second-order transitions, whereas a phenomenological Lan- 
dau mean-field theory predicted a richer structure. 19 The Lan- 
dau theory also suggested a possibility of a single second- 
order transition. 19 We simulated the system for J z = 1, 0.25, 
0.1, 0.0625, 0.04, 0.0278, and 0.0204. The aspect ratios of 
the lattices are r = 1, 1/2, 1/3, 1/4, 1/5, 1/6, and 1/7, re- 
spectively, and were determined to reduce finite- size effects. 

In Fig. [3 we show the internal energy distribution at a tem- 
perature around which the Binder parameters show character- 
istic features of a phase transition (a shift from the high-T 
Gaussian value to the trivial low-T value) for several values 
of J z . The distribution exhibits bimodal structure, which is an 
unambiguous signature of a first-order transition. The peak- 
to-peak distance, AE, is a finite-size estimate of the latent 
heat. As shown in the inset of Fig. [3 AE increases monoton- 
ically as a function of J z in this region. 

It is natural to ask what are the broken symmetries be- 
low the first-order transition. In Figs. H3a) and [TUlb) we 
show the temperature dependence of the Binder parameters 
U a and U m for J z = 0.0625 and 0.0278, respectively. Both 
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T T 

FIG. 1 1 . (Color online) Two possible scenarios that may explain our 
observations in the intervening region are shown, (a) The two transi- 
tions continue to be of second order until they collapse into the direct 
first-order transition, (b) Only the lower-T 0(3) transition becomes 
of first order before the transitions merge. The arrows schematically 
indicate the possible ways in which the system is cooled down for 
0-01 %Jz< 0.0204. 



of them show diverging behavior, indicating discontinuity of 
the corresponding order parameters 3 The same behavior is 
observed for the J z values listed above, except for the case 
J z = 0.0204, which will be discussed later. The peaks 
are sharper for larger J z values. Based on this observation 
and the monotonically increasing value of AE, we conclude 
that the 3D Ising-0(3) model with a large enough inter-layer 
coupling undergoes a single first-order transition. Our con- 
clusion agrees with the RG analysis of i^eff* 36 while it dis- 
cards the other previous scenarios, namely the single second- 
order transition suggested by the Landau mean-field theory 19 
or the always separate second-order transitions suggested by 
the previous large-TV approximation. 10 However, recently an- 
other large-TV approach yielded a phase diagram showing a 
single first-order transition when the Ising and the 0(3) tran- 
sitions are merged, in agreement with our results, 43 On the 
other hand, the failure of the Landau mean-field theory at this 
point is not surprising because the system is below the upper 
critical dimension d = 4. 

Finally, we briefly discuss the region where the two tran- 
sitions merge. As we mentioned before, there is some level 
of uncertainty in this intervening region. By comparing 
Fig. 03 a) andfTOtb), we notice that there is a monotonic ten- 
dency in the peak structure. As J z decreases, the peak of U a 
is drastically suppressed as compared to that of U m . Indeed, 
as shown in Fig. H3c), U a for J z = 0.0204 does not ex- 
hibit an evident diverging behavior for the explored system 
sizes, while U m clearly does. One possible explanation is that 
finite- size effects smear the discontinuity of the Ising order 
parameter because the Ising correlation length is larger than 
L (not shown). Thus, based on our numerical results, we can- 
not discard a scenario in which the direct first-order transi- 
tion to the lowest-T ordered phase ends up at a critical point 
where it starts splitting into two second-order transitions [see 
Fig- Eta)]. However, our RG analysis indicates that such a 
critical end point would not be stable. The most plausible 
scenario corresponds to the existence of a finite region where 
the system in the lowest-T phase first recovers the 0(3) sym- 
metry via a first-order transition, while the symmetry is 



recovered at a higher temperature via a second-order Ising- 
like transition [see Fig.QjJb)]. This implies that U m should 
diverge around the first-order transition, whereas U a should 
reach a trivial and finite low-T value below T = T c \ for large 
enough system sizes. Although we still do not have enough 
evidence to confirm this scenario, our results suggest that it 
may occur near J z = 0.0204. 



V. SUMMARY 

In summary, we have studied the Ising-0(3) model on 
quasi-2D lattices. This is an effective model of the J1-J2 
Heisenberg model in which J2 is dominant. By solving 
this effective Hamiltonian, we identified the region where the 
0(3)- symmetric Ising ordered phase exists. For sufficiently 
large inter-layer coupling, we found that a single first-order 
transition occurs between the paramagnetic phase and the 
lowest-T ordered phase, in agreement with the previous RG 
treatment on H t a |2S Although the question of how this first- 
order transition splits into two transitions remains as an open 
problem, the scenario shown in Fig. fTTT b) provides the most 
reasonable explanation of our numerical results. 

Our results provide a qualitative explanation for the se- 
quence of transitions observed in ferropnictides as a func- 
tion of the ratio between the inter-layer and intra-layer ex- 
change couplings. According to these results, the separate 
structural and SDW transitions observed in the quasi-2D 1111 
compounds 4 ^ are caused by the fragility of the continuous 
SDW order against fluctuations, which makes it more sensi- 
tive to the magnitude of the inter-layer coupling. The sug- 
gested structural transition in Li2VOSi04— is also naturally 
explained by the same mechanism. The first-order nature of 
the simultaneous structural and SDW transition observed in 
most of the more 3D 122 parent compounds^S is also con- 
sistent with our results. This is related to the absence of the 
stable RG fixed point^ Remarkably, a recent measurement 
also found that a sequence of transitions that is entirely con- 
sistent with the scenario presented in Fig. fTTT b) takes place 
in BaFe2As2*22i2! Therefore, despite the oversimplified nature 
of our local-moment model for a microscopic description of 
the metallic ferropnictides, we have reproduced their qualita- 
tive phase diagram (Fig.0. Our results thus indicate that the 
J\- J2 Heisenberg model and the related Ising-0(3) model are 
good starting points for describing the universal properties of 
these compounds. 
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Appendix A: EFFECTIVE HAMILTONIAN FOR THE 
ISING-O(iV) MODEL 

Let us consider a generalization of the Ising-0(3) model 
replacing the three-component spin by the O(N) spin: H = 
— Jij (1 + <J i (J j) S« • Sj, where oi and S$ are the Ising 
and the O(N) spins, respectively. In the following, we derive 
an effective Hamiltonian for this model. By using a Gaussian 
transformation, we introduce the auxiliary A^-component vec- 
tor fields 4>a ~ S and 4>b ^ crS. By introducing = (3Jij, 
we obtain: 



e~ m = exp Q Si ■ KijSj^J exp Q cr^ • K^ajS, 
ex J D[(f) A }D[(f)B} exp ' K ij 1( l>A,j + • S, 

x exp \-\(t>B y i • K^^bj + <f> B ,i • cTiSi 



exp[(0 A ,i +CT0B,i) • 



(Al) 



where the summation rule for duplicate indices is assumed. 
Since we have decoupled the spins on different sites, we can 
trace them out on each site: 



Z oc D 



I 



[(j) A ]D[(j) B } exp | -- 0<M • K ij 1( t>a,j J 

\ a—A,B J 

x Y[ Tr s,a exp [(4>A,i + CT0B,i) • S] 
D[4> A ]D[4> B ]expl~ Yl ^i-Kr/^j] 

\ a—A,B J 



exp < Y ln 



a=±l 



(A2) 



Here, 



Giv(j)=Trs exp (j • S) 

J^^(5(S 2 -l)exp (j-S) 
" f d N S 6 (S 2 -l) 

f -xn r(AT/2) 



n=0 



2 2n n j r(iV/2 + ra) 



(A3) 



is a single- site generating function. The following terms ap- 
pear in the expansion of y }2 a=±1 Gn ((/>a, % + &(f>B,i)'- 



C2 (\4>A,i + 0B,i| 2 + \4>A,i ~ <t>B,i\ 2 ^) 

= 2c 2 (|0^| 2 + |0b^| 2 ) 
= 2c 2 $ 2 , 



(A4) 



(|0 a ,,| 2 + |0b,,| 2 ) 2 +4(0 a ,,.0 b ^) 2 



2c 4 



2c 4 $ 4 . 



(A5) 



By truncating at fourth order in powers of the fields 4> A and 
<f>B, we obtain: 



ln 



^ G N ((j) A ,i + CT0B,i) 



=±i 



ln(ci+2c 2 $2 + 2c 4 <l> 4 + ...) 

ci ci 2 V ci ' 



(A6) 



From this expression we obtain 

Z = J D[c/> a }DIcI>b} exp (S[4> A , 4b]) , (A7) 

with 

S[</>A, 0b] 

=i x: *-,*-^v«j-E^(i^.«i a+ i^.«i a ) 

+ ^ U (\4> A ,i\ 2 + |0B,i| 2 ) + A(0A,i * 

(A8) 

Here, one can verify that the coefficient u = 2 (C2/C1) 2 — 
2c 4 /ci is positive and A = — 8c 4 /ci is negative. It is straight- 
forward to rewrite the quadratic terms in the form given in 
Eq. ©. 



Appendix B: REVIEW OF THE RG TREATMENT ON THE 
EFFECTIVE HAMILTONIAN 



In this appendix, we derive the one-loop RG flow equa- 
tions of H Q ff [Eq. [3) . The final result was first presented by 
Aharony^ in the rather different context of amorphous mag- 
nets. For completeness, we will consider the generalization to 
O(N) of our 0(3) invariant Hamiltonian H Q a (N is the num- 
ber of components of each spin). 

We start by analyzing the stability of the so called decou- 
pled fixed point (DFP). This is the Wilson-Fisher fixed point 
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with O(N) symmetry and u* = e/ [8 (TV + 8)] + O (e 2 



1 AB 



A* 



0. The reason we are interested in this fixed 



point is two-fold. In the first place, the values of the bare cou- 
pling A or uab can be very small for some frustrated magnets 
such as the J1-J2 model with J 2 ^> J\ or the quasi-2D bet 
lattice model.^^"— In the second place, we can discuss the 
stability of the DFP very accurately using a non-perturbative 
scaling argumen t 25 ! 36 ! 47 because the sublattices are trivially 
decoupled and the 0(7V) fixed point has been studied very 
extensively (there are very accurate estimations of the corre- 
sponding exponents). 

The stability of our DFP is determined by the scaling di- 
mensions of uab and A, which can be obtained from the two- 
point correlators of the conjugate scaling operators: First, 

<|<M*l)| 2 |<M*l)| 2 \<f>A(X2)\ 2 \(f>B(X2)\ 2 ) D 



(\ < j )A (x 1 )\ 2 \ ( j )A (x 2 )\ 2 ) D (\cl )B (x 1 )\MB(x 2 )\ 2 ) 1 



OC \X\ - x 2 \ 



-Ax t 



(Bl) 



Here the average {•) D is taken under the condition uab = 
A = 0, and x t = d — \ jv is the scaling dimension of the 
energy-density operator at the 3D-0(iV) DFP. Equation (IB 11) 
shows that the scaling dimension of | 4>a | 2 1 4>b | 2 is simply 2x t 
and thus the RG eigenvalue is 



t>UAB 



d — 2x+ 



2 -dv 



a 



(B2) 



where we have used the hyperscaling relation a = 2 — dv. 
Because the specific-heat exponent a of the 3D-0(7V) models 
is known to be negative for TV > 2, we conclude that y^F AB < 
0, i.e., the DFP is stable against the uab term. 

Now we discuss the relevance of the A term. By introducing 
the traceless symmetric quadrupolar tensors Q% v = (j)%(j) v a — 
N~ 1 S^ u \4> a \ 2 (a = A, B), we can decompose the A term in 
the following way: 

(<Pa ■ 4>b? = Q7Q7 + ^l<^| 2 |<H 2 - (B3) 

By using the 0(7V) invariance of the decoupled Hamiltonian 
we obtain: 



(<f>A • 4>b) 2 Oi) {<f>A • 4>b) 2 (x 2 ) 



D 



= (Q7(xi)Q A x (x 2 )) D (Q'b(xi)Qb k (x 2 )) d 
+ 7V- 2 (|^( a;i )| 2 |^(x 2 )| 2 ) I) <|0 B ( a;i )| 2 |0 i3 ( a;2 )| 2 ) £) 
Cqq Cu 



xi - x 2 \ Axq \x\ - x 2 



\4x t ' 



(B4) 



Here Cqq and C u are nonzero coefficients and xq is the scal- 
ing dimension of the quadrupolar order parameter. Since the 
second term is irrelevant at the DFP, the scaling dimension of 
the A term is equal to 2xq. By defining vq = d — xq, we 
obtain 



y{ D) =d-2x Q = 2y Q - d. 



(B5) 



Reference 48 provides estimates of vq (denoted as y 2 there) 
for TV = 2, 3, 4, 5 and 16. The value of vq = 2 for N — >> 00 is 



also provided. In all of these cases we find y[ D) > 0, meaning 
that the DFP is unstable in the presence of the A term. 

It is then natural to ask whether a stable fixed point exists in 
the proximity of the unstable DFP. In the following, we show 
the results obtained by expanding around the Gaussian fixed 
point in 4 — e dimensions to O(e). Here, we use the notation 
introduced by Cardy 49 and derive the flow equations to O(e) 
by applying the operator-product expansion (OPE) method. 
We assume that the operators are normalized in such a way 
that (0£ (xi) <\) v h (x 2 )) = Subtil \xi ~ x 2 \~ {d ~ 2) at the Gaus- 
sian fixed point (a, b = A, B and 1 < /i, v < N) and that 
they are normal-ordered in a sense described in Ref . |49|. The 
following OPE's are sufficient to construct the RG equations: 

Vv • Vv = 47V + 4Vv + + 2^ nAS , 
ipr'ipu =4(7V + 2)Vv + 8^n, 
■ ^u AB = 47V^r + Sif) UAB j 

^r'^X= 2?/V + 8^A, 

^u'^u = 247V 2 + 32 (TV + 2) + 8 (TV + 8) i/> u , 
^w^u AB = 8(7V + 2)^ AB , 

^u'^X= Hu AB + 16^A, 

^u AB ■ = 4iV 2 + 87V Vv + 2Ni(j u + 16^n AS , 
*p UAB ■ = 47V + 8Vv + 2^n + 16^ A , 

V^a • = 47V 2 + 4 (TV + 1) t/v + 2^n + ^u AB 

+ 4(7V + 2)^A- (B6) 

Here, Vv = \cf) A \ 2 + \(f> B \\ ^ = |0a| 4 + \^b\\ ^ Uab = 
|0a| 2 I <\>b 1 2 and = ((f) a -0b) 2 are short-hand notations 
for the scaling operators. The RG flow equations to O(e) are 
entirely determined by these OPE coefficients* 4 ^ 

dv 

— = 2r - S(N + 2)ru - 4rA - 4Nru A B - • • • , (B7) 
du 



^■ = eu-8(N + S)u 2 - 2A 2 - 4u AB X 



2Nu 



AB 



(B8) 



du 



AB 



dl 
dX 



eu A B ~ 16?xA - 16(iV + 2)uu A b ~ 4A 2 



16^ B - . . . , (B9) 



^ = eA - 32u\ - 4(7V + 2)A 2 - 32u AB \ 

(BIO) 

where we assume that the fixed points of physical interest are 
located in the region where r = 0(e 2 ), u = O(e), uab = 
0(e) and A = 0(e). 

Let us first discuss the physically relevant cases N = 2,3. 
The fixed points to O(e) for N = 2 are as follows: 

• Gaussian fixed point: (u,uab,X) = (0,0,0) 

• XY DFP: (u,uab,X) = (e/80,0,0) 

• 0(4)-like: (u,u A bA) = (e/96, e/48, 0) 
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• (u,u AB ,X) = (e/160,3e/80,-e/40) 

• (u,u A bA) = (e/160,e/80,e/40) 

The fixed points to O(e) for TV = 3 are as follows: 

• Gaussian fixed point: (u, uab, X) = (0, 0, 0) 

• 0(3) DFP: (u,u A b,X) = (e/88,0,0) 

• 0(6): (u,uab,X) = (6/112,6/56,0) 

• (u, uab, A) = (3e/272,e/136,0) 

• (u,u A b,X) = (e/136,e/68,e/68) 

• (u,u AB ,X) = (e/176,e/88,e/44) 

The most important conclusion is that none of these fixed 
points is stable to O(e). Therefore, this simple RG calcula- 
tion suggests that the biquadratic coupling A between the two 



O(N) subsystems leads to a first-order phase transition^ It 
is interesting to note that, in contrast to the result obtained 
by directly evaluating the correlation function, the one-loop 
expansion indicates that uab is a relevant perturbation at the 
DFP for TV < 4 and d < 4: 

V&l = ( 4 " N ) e l ( N + 8 ) + °^ 2 )' (Bll) 

Note that y^ B is positive to 0(e) for TV < 4 and d < 4, while 
our Eq. flSJ shows that < for TV > 2 in d = 3. This 
discrepancy must be eliminated by the higher-order terms of 
the e expansion. Although the concomitant change that will 
appear in the RG flow structure is unclear, numerical studies 
of microscopic Hamiltonians, such as the Ising-0(3) model in 
the present work or the coupled XY model in Ref. l27l confirm 
the absence of a stable fixed point. The model does not have 
a stable fixed point in the region A < even for larger values 
of 
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